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A  Practical  Approach  to  Karmarkar’s  Algorithm 
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Department  of  Operations  Research 
Stanford  University 
Stanford.  CA  94305 


ABSTRACT 

A  practical  approach  to  implementing  Karmarkar's  algorithm  is  discussed. 
A  variant  of  the  algorithm  is  proposed  which  still  has  polynomial  complexity 
and  which  eliminates  the  need  for  Karmarkar's  canonical  form.  This  method 
allows  upper  and  lower  bounds  to  be  used  and  does  not  require  knowledge  of 
the  objective  value.  Some  heuristics  are  given  which  alleviate  certain 
computational  difficulties  that  arise  when  a  practical  implementation  of  the 
algorithm  is  attempted.  A  FORTRAN  program  is  described  that  allows  one  to 
study  its  convergence  properties. 
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Narendra  Karmarkar  of  AT&T  Bell  Laboratories  has  described  a  new 
polynomial-time  algorithm  for  linear  programming  that  has  generated  a  large 
amount  of  press  coverage.  The  important  theoretical  result  he  presented  at 
the  Symposium  on  Theory  of  Computing  in  Spring  of  1984  has  been  well 
received,  but  statements  regarding  the  practical  efficiency  of  the  method  have 
stirred  up  much  controversy  among  experts  in  the  field.  His  initial  claim  that 
the  algorithm  was  50-100  times  more  efficient  than  the  simplex  method  was 
recently  repeated  by  him  at  the  Operations  Research  Society  of  America 
meeting  in  Boston  in  April,  1985,  but  the  type  of  problems  for  which  these 
results  were  claimed  were  for  a  very  small  number  of  test  problems  specially 
structured  to  favor  his  approach  151.  Tomlin  has  attempted  to  solve  some  less 
structured  problems,  and  has  not  been  able  to  duplicate  Karmarkar’s  less 
spectacular  computation  times  on  these  smaller  problems  [81.  Karmarkar  [41 
has  claimed  that  ‘any  undergraduate  should  be  able  to  read  my  paper  and  write 
a  program  that  is  10  times  faster  than  the  simplex  method.*  This  manuscript 
describes  the  efforts  of  the  author,  a  graduate  student,  to  validate  this  last 
claim. 

Section  1  describes  one  possible  way  to  apply  Karmarkar 's  method  to 
general  linear  programs  and  some  of  the  basic  difficulties  that  arise.  Section 
2  proposes  a  variant  of  the  method  that  remains  polynomial  and  allows  general 
linear  programs  to  be  solved.  Section  3  discusses  how  the  method  is  applied 
in  practice,  the  generation  of  interior  feasible  pofnts,  and  the  problem  of  null 
variables.  Section  4  describes  an  iterative  scheme  to  find  the  projected 
gradient  on  each  iteration.  Section  5  discusses  an  implementation  of  the 
algorithm  and  some  computational  results.  Section  6  concludes  the  paper  with 
a  discussion  of  some  future  research. 
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An  initial  reading  of  Karmarkar's  paper  [2]  indicates  that  a  practical 
Implementation  of  the  algorithm  is  not  easy  to  come  Dy.  One  difficulty  arises 
because  projective  transformations  are  used  to  transform  any  general  linear 
program  to  the  following  canonical  form  with  a  homogeneous  right-hand  side 
and  a  convexity  constraint: 

min  c^x  c  € 

(LPO)  s.t.  Ax  =  0  A  £  b  £ 

ejx=  1  e,(  =  (t.l . 1)'^£R'^  .  Vk>  1 

X  >  0. 

Todd  and  Burrell  17]  and  Tomlin  [8]  have  described  methods  to  do  this 
transformation,  but  these  methods  either  add  unnecessary  rows  and  columns  to 
the  linear  program,  and/or  depend  on  the  existence  of  some  upper  bound  on  the 
sum  of  all  of  the  variables.  For  practical  implementations,  finding  such  a 
bound  may  be  as  difficult  as  solving  the  linear  program  and  using  this  bound 
can  cause  numerical  instability  during  the  course  of  solution.  The  procedure 
described  herein  has  the  advantages  that  only  one  dense  column  in  Phase  I 
needs  to  be  added  and  the  bad  scaling  properties  of  previously  suggested 
methods  are  not  present. 

Suppose  we  are  given  the  following  linear  program: 

min  c^x  c  c  R'^ 

(LPl)  s.t.  Ax  =  b  A  f  b  f  R"' 

X  >0 

with  known  feasible  point  w  £  R^,  w>0  and  known  minimum  objective  value  z**. 

Karmarkar  [4]  Indirectly  suggested  that  a  projective  tranformation  be  used  to 
bring  this  program  Into  the  canonical  form  of  (LPO).  thereby  creating  the 
following  linear  program: 


(LP2) 


min  [c'D.-z^lir 
s.t.  [AO.-b]>r  =  0 

=  I 

>r>  0. 


where  D  =  diag[w^.  W2 . w^j.  (LP2)  is  in  the  canonical  form  for 

Karmarkar's  algorithm  (2)  (henceforth  KA),  It  also  has  not  changed  the 
sparsity  structure  of  the  original  problem,  except  for  the  addition  of  the 
convexity  constraint.  If  the  KA  algorithm  is  applied  directly  to  this  linear 
program,  a  new  projected  space  is  created  on  each  iteration  in  order  to  move 
the  current  iterate  to  the  center  of  a  simplex.  Therefore,  a  series  of  three 
linear  spaces  are  used  to  solve  the  linear  program^  The  first  space  is  that  of 
(LPl),  the  second  is  (LP2).  and  the  third  is  artificially  created  on  each 
iteration  when  the  projective  step  is  made.  Todd  and  Burrell  [7]  and  Tomlin  [81 
have  not  used  the  form  of  (LP2)  to  solve  (LPl),  and  the  equivalence  of  (LPl) 
and  (LP2)  is  not  apparent  in  Karmarkar's  paper  [21.  due  to  the  addition  of  the 
column  for  "Xp^+j. 

Iterates  of  the  KA  algorithm  are  calculated  in  terms  of  "x,  a  vector  In 
When  the  algorithm  stops  with  a  solution  1?*,  the  solution  x**  to  (LPl)  is 
given  by  the  inverse  projective  transformation  x*'=[D,0l>^/j^^  I .  Furthermore, 

the  knowledge  of  the  finite  minimum  to  (LPl)  guarantees  the  existence  of  >?*. 
with  an  objective  value  for  (LP2)  of  0,  as  the  following  theorem  shows: 

Theorem  1  Suppose  (LPl)  has  a  finite  optimal  solution  x”.  Then  (LP2)  has  an 
optimal  solution  at  >?*  =  (D"’x*',t)/(1+ej!jD'’x*‘). 


Proof  It  is  clear  that  x**  Is  feasible  for  (LP2)  and  that  Its  objective  value  Is 
0,  since  c^x**  =  z**.  Now  suppose  3  ^  that  is  feasible  for  (LP2)  and  that 
[c^D.-z**!^  <  0.  There  are  2  cases  to  consider,  depending  on  the  value  of  j . 

Case  1:  Assume  Xp*|>  0.  Let  x=[D,0]x/ Xp+ 1 .  Then  x  is  feasible  for 
(LP2).  [c^D,-z**lx  <  0  =>  c^[D,0]x  <  z^Xp+j  ^  c^x  <  z"  =>  x"  is  not  optimal  for 


Case  2=  Assume  1  =  0.  Lei  x=x*‘+lD.0]x.  Then  [AD,-blx=0  =>  AlD,0lx=0, 
so  that  X  is  feasible  for  (LP1).  Ic^D.-z^lx  <  0  =>  c^tD.Olx  <z*‘Xp+|  => 

c^(D.0]x<0  =>  c^x=z*‘*c^lD,Olx  <  z**  =>  x**  is  not  optimal  for  (LP1). 

Hence,  in  either  case.  cannot  exist,  and  >?*  must  be  optimal  for  (LP2).  • 

In  fact,  one  can  see  from  the  above  proof  that  x^+pO  and  [c^D.-z’^lx  <0 

implies  the  existence  of  an  unbounded  optimal  solution  in  (LP1)  by  simply 
moving  in  the  direction  ID.Olx  from  x^  In  the  feasible  region  for  (LPI). 

Difficulties  with  running  the  KA  algorithm  directly  on  (LP2)  are  due  to 
two  factors.  If  the  optimal  >?*  found  by  the  algorithm  has  =  0.  then  a 

nonnegative  objective  value  for  x**  only  Implies  that  the  original  problem  has 
an  infinite  ray  along  which  It  is  optimal.  Unfortunately,  the  KA  algorithm  does 
not  yield  the  endpoint  of  this  ray  for  (LPI).  It  seems  difficult  to  prove  that 
this  condition  will  not  occur.  The  second  difficulty  arises  when  one  tries  to 
extend  the  basic  idea  of  projecting  (LPI)  to  (LP2)  when  the  linear  program  has 
variables  with  upper  bounds.  This  problem  xcurs  because  the  upper  bounds  do 
not  project  to  a  simple  linear  ratio  test  when  performing  the  iterative  step. 
The  projective  transformation  used  on  each  iteration  that  maps  (LP2)  to  the 
third  space  mentioned  above  destroys  this  linearity.  With  this  motivation,  we 
modify  the  basic  idea  of  Karmarkar’s  algorithm  to  avoid  these  and  other 
difficulties. 

Intuitively,  projecting  the  current  iterate  to  the  center  of  the  feasible 
region  is  also  a  good  way  of  solving  linear  programs.  In  fact,  there  seems 
good  reason  to  consider  an  algorithm  which  simply  brings  each  current  iterate 
into  the  center  of  some  simplex  by  a  projective  transformation,  moves  in  the 
direction  of  optimality  of  the  linear  function  over  the  inscribed  sphere,  and 
then  uses  an  inverse  projective  transformation  to  bring  the  new  point  back  to  a 
point  in  the  original  space.  Such  an  algorithm  works  only  with  iterates  in  the 
original  space  of  the  linear  program  and  uses  projective  transformations  only 
to  move  from  the  current  iterate  to  the  next  one.  When  doing  this,  there  is  no 
need  to  bring  a  linear  program  to  any  canonical  form  like  (LPO).  Furthermore, 
there  are  now  only  two  spaces  used:  the  original  linear  program  (LPI)  and  the 
space  created  on  each  Iteration.  The  next  section  will  show  that  this  variant 
of  the  KA  algorithm  also  remains  polynomial. 


For  this  algorithm,  we  assume  we  are  given  a  linear  program  in  the  form 
of  (LPl),  the  feasible  point  w.  and  the  known  objective  value  z".  Later,  we 
will  indicate  how  such  a  point  w  can  be  obtained  and  how  z"  need  not  be  given. 
The  NPA  algorithm  is  as  follows: 

begin 

X^«-W  ; 

k  ♦-  1  : 

while  (not  optimal)  jlfi 

hfigio 

D  «-  diagtx'j.x^ . xj^.H 


fiOSt 


(A.-b)D 


Cf.*-  [c.-z**] 

Cp^  (l-B'^(BB'’’)"'BlC)Cr 
c  ^  Cp/||Cp| 

b  <-  (Db')/(eJ,|Db‘) 

B  b/b^^, 

k  k+l 

x^*-  b|j  (i.e,  the  first  n  components  of  b) 

end 
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For  theoretical  purposes,  the  NPA  algorithm  will  be  considered  terminated 
at  optimality  when  (c^X“Z*‘)<2'^.  where  q=0(L).  as  in  Karmarkar's  paper.  Also, 
r  is  set  equal  to  l/yn(n*l).  which  is  the  radius  of  the  largest  inscribed  sphere 
in  a  simplex  of  dimension  n.  The  parameter  «  is  set  to  .25.  to  facilitate  a 
polynomial  proof  of  convergence.  In  a  following  section,  it  is  shown  how  the 
above  algorithm  is  modified  for  practical  purposes  by  selecting  <x  differently 
and  using  an  appropriate  convergence  criterion. 

In  order  to  show  the  algorithm  is  still  polynomial,  it  is  necessary  to 
invoke  the  use  of  a  potential  function  similar  to  the  one  used  by  Karmarkar.  In 
fact,  this  potential  function  will  relate  in  a  special  way  to  that  of  Karmarkar. 
We  define  the  potential  function  as: 

n 

g(x'^;c;z*‘)  =  (n+l)itn(c^x’^-z*‘)-  2  Mx^)  • 

The  following  lemmas  will  show  how  this  potential  function  leads  to  the 
polynomial  complexity  of  the  NPA  algorithm. 

Lemma  1  g(x*^*  ^c;z**)<  g(x*^;c;z'*H,  where  8=1/8,  as  in  131. 

Proof  In  each  Iteration  of  the  NPA  algorithm,  the  current  iterate  x*^  is 
projected  to  the  center  of  a  simplex  to  create  a  linear  program  in  the 
canonical  form  of  (LPO)  with  n*l  variables,  and  a  known  minimum  of  0,  by 
Theorem  1.  This  linear  program  has  the  form: 

c^=lc^.-z**lD 
A=(A.-blD 


min  £^x’ 
(LP3)  s.t.  Ax'=0 


x’>0 


Applying  Theorem  2  of  (3],  one  finds  that  f'(b')lf’(ao)-5,  where  b’  is  as  in 
the  NPA  algorithm,  ao=ep,*l,  8=1/8,  and  f  is  Karmarkar’s  potential  function  for 

the  transformed  linear  program  (LP3),  as  found  in  [31: 


7  TJ  7  ? 


r^  *  2  W'j  ' 


n+1 

f*(x')=(n+l)iln(c^x*)-  2  • 

j=1 

Hence,  the  reduction  in  Karmarkar's  potential  function  implies  that 

n+1  n+l 

(n+OMc^b')-  2  Mb j)  <  (n*  DMc^ao)-  2  Anj^  -  8  . 
j=1  j=l 

Since  B  =  (Db')/bJ^+ j.  effectively  scaling  B  so  that  it  is  feasible,  we  can  write 

n+l 

(n+  OMc^B-z**)-  2  Mbj/Dj  j)  <  (n+  l)MDc^epj-2*‘)-8 . 
j=l 

Now,  Br,4)=l  and  D=diaglx*5.x|.....xj^,ll,  and  x‘^*'=B({ . together  imply 

n  n 

(n+ 1  )Mc^x’^^ ' -z")-  2  Mxj* ' )  <  (n+ 1  )Mc^x’^-z*‘)-  2  Mx*!)  -  5  . 

j=l  j=l 

which  yields  the  final  result,  by  the  definition  of  g(x*^:c;z*’).  • 

lemma  2  g(x'^:c;z*‘)<  g(w;c:z*‘)-k8. 

Proof  Obvious,  by  Lemma  1 .  • 

Lemma  3  After  k  iterations  of  the  modified  algorithm. 

(c^x'^-z**)  <  exp(q-  ^Kc^w-z^Xdet  diag[w^,W2.....Wpjl)^*'^^^‘^ . 

Proof  Lemmas  2  and  1  together  imply 
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n  n 

(n+OMc^x^^-z**)-  2  Mx*-)  <  (n*1)jln(c^w-z’*)- 2  MWj)  -  k8  . 

j=1  j=l 

Since  we  assume  there  to  be  an  optimal  solution  to  the  linear  program,  each 
component  of  x*^  is  bounded  above  by  e^.  Hence, 

n 

Mc^x'^-z**)  <  +  AnCc^w-z*)  -  ^  Mwj)  -  ^  . 

After  exponentiating  both  sides,  the  desired  result  is  obtained.  • 

One  may  argue  that  the  exponential  term  in  the  result  of  Lemma  3  is  quite 
large.  However,  as  k  increases,  this  term  will  eventually  become  small,  and 
this  yields  a  proof  that  a  polynomial  number  of  Iterations  will  occur. 
Furthermore,  the  other  two  terms  in  the  product  are  of  the  same  order  as  the 
exponential  term  (i.e.  their  logarithms  are  polynomial  in  the  size  of  the  input) 
and  these  terms  do  not  effect  the  overall  complexity  of  the  NPA  algorithm. 
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The  first  major  modification  to  the  NPA  algorithm  is  in  the  choice  of  «. 
Karmarkar  [4]  had  originally  suggested  that  the  following  ratio  test  be 
performed: 

(x^  -  max  {l/(n+1)rcj  .  Cj  >  0  } 

This  would  find  a  blocking  variable  Xj.  (We  know  that  oc**>]  because  (x‘‘=l 

corresponds  to  a  distance  of  one  radius  of  the  inscribed  sphere.)  One  would 
then  use  o<=po<**  as  the  step  length  in  the  algorithm,  where  p=.85.  .90.  or  .95. 

say.  Then  p  represents  a  fraction  of  how  close  one  would  get  to  any  face  of 

the  positive  orthant  in  R^.  Tomlin's  initial  experimentation  indicated  that  this 
was  a  good  approach,  as  values  of  <x**>l  were  often  observed,  which  indicated 
that  the  algorithm  was  moving  in  profitable  directions  to  points  far  outside 

the  sphere  inscribed  in  the  simplex.  The  cost  of  computing  the  value  is 

simply  n  ratio  tests.  Todd  and  Burrell  [71  have  suggested  that  a  line  search  be 
performed  in  the  direction  c  to  minimize  the  potential  function.  It  is  not 
clear  how  beneficial  this  search  would  be,  versus  the  above  heuristic,  which 
seems  to  work  well  in  practice. 

Given  a  general  linear  program  in  the  form  of  (LPl),  the  first  problem  that 
arises  is  that  of  finding  an  interior  point.  As  suggested  by  Karmarkar  in  [21, 
one  may  set  up  the  following  linear  program  to  get  an  interior  point: 

min  X 

(LP4)  s.t.  Ax+(b-Aep)X  =  b 
x>0.  X>0  . 


This  linear  program  has  an  optimal  value  of  0  if  and  only  if  (LPl)  is 
feasible.  Furthermore,  the  vector  x=erj.  X=l,  is  feasible  for  this  linear 

program.  We  can  then  apply  the  NPA  algorithm  directly  to  this  linear  program 
to  yield  a  feasible  point.  If  the  program  (LPl)  is  infeasible,  we  will  probably 
take  many  iterations  to  detect  that  condition.  However,  if  (LPl)  has  a  strictly 
interior  point,  X  will  become  a  blocking  variable  in  relatively  few  iterations, 
with  the  Iteration  count  depending  upon  how  far  e^  lies  from  any  solution  of 

(LP4).  When  X  blocks  in  the  ratio  test  mentioned  above.  «  is  set  to  oc**,  and 
the  algorithm  obtains  a  strictly  interior  point. 
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A  heuristic  argument  can  be  given  as  to  why  the  above  procedure  yields  an 
interior  point,  if  one  exists.  In  each  iteration  of  the  algorithm,  the  variable  X 
is  projected  to  a  simplex  as  one  of  the  variables  of  the  linear  program.  Since 
X=0  corresponds  to  all  feasible  solutions  of  (LPl)  and  to  all  optimal  solutions 
of  (LP4)  and  also  to  a  face  of  the  simplex  artificially  aeated  in  each  iteration, 
one  can  see  that  the  linear  programs  (LP3)  and  (LP4)  have  multiple  optimal 
solutions.  The  KA  and  NPA  algorithms  both  tend  to  go  to  the  "center*  of  such 
a  face  of  multiple  optima,  when  starting  with  a  point  sufficiently  far  away 
from  any  boundary.  This  Is  a  result  of  minimizing  over  the  sphere  inscribed 
Inside  the  simplex,  and  the  fact  that  this  sphere  Is  tangent  to  the  face  X=0  of 
that  simplex. 

In  Phase  I.  problems  arise  In  the  presence  of  null  variables.  For  (LPl),  Xj 
Is  said  to  be  a  null  variable  If  Xj=0  in  all  feasible  solutions  for  (LPl).  Hence, 

a  linear  program  with  any  null  variables  has  no  strict  interior.  It  can  be 
shown  that  a  linear  program  has  one  or  more  null  variables  if  and  only  if  every 
basic  feasible  solution  of  (LPl)  Is  degenerate.  One  should  also  note  that 
randomly  generated  linear  programs  rarely  have  null  variables,  but  null 
variables  seem  to  occur  quite  often  in  practical  problems. 

The  behavior  of  both  the  KA  and  NPA  algorithms  on  (LP4)  is  quite 
interesting.  The  solution  of  (LP4)  still  lies  on  the  face  X=0,  but  this  face  is 
not  n-dimensional,  as  in  the  case  when  no  null  variables  are  present.  It  has 
been  observed  on  various  test  problems  that  X  does  not  become  the  blocking 
variable  In  a  small  number  of  iterations;  Instead,  the  algorithm  tends  to  "tall 
off"  as  it  converges  to  a  solution  with  X=0.  However,  as  Iterations  progress 
the  algorithm  seems  to  identify  the  null  variables  in  an  Interesting  way.  Since 
the  algorithm  Is  converging,  the  change  In  those  variables  that  are  not  null  is 
very  small  compared  to  those  that  are  null.  If  x*|  is  small,  it  must  be 

projected  relatively  far  to  the  center  of  the  next  simplex  in  the  next  iteration. 
It  has  been  observed  that  for  null  variables  that  are  converging  to  0,  the  value 
of  Xj  Is  very  small  compared  to  that  of  its  projected  gradient  3j. 

A  heuristic  which  takes  advantage  of  this  observation  has  been 
implemented.  If  a  variable  Xj  is  the  blocking  variable  in  the  ratio  test  and 

(xj/(cj/Xj))  <  e,  then  «  is  set  to  «•*,  driving  Xj  to  zero.  Note  that  the  test 

divides  Cj  by  Xj  so  that  the  values  compared  are  related  to  each  other  within 

the  space  of  (LPl).  For  practical  testing,  a  value  of  was  found  to  be 
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suitable.  After  that  variable  has  been  set  to  0,  it  is  effectively  ignored,  and 
the  algorithm  prxeeds  until  X  becomes  blocking  in  a  ratio  test.  During  that 
time,  other  null  variables  may  be  identified  and  set  to  0  as  well.  As  soon  as 
an  interior  point  (ignoring  any  null  variables  that  have  been  identified)  has 
been  generated,  one  can  then  solve  (LP1). 

The  final  difficulty  in  solving  the  linear  program  is  the  requirement  that 
z**.  the  optimal  objective  value  of  (LPl).  be  input  to  the  algorithm.  It  does  not 
seem  practical  to  solve  the  combined  primal-dual  program  as  suggested  by 
Karmarkar  in  [31.  Yinyu  Ye  has  used  in  [91  a  "cutting  objective"  method  to 
solve*  this  problem.  This  method  is  extremely  easy  to  apply  to  the  NPA 
algorithm,  as  it  simply  sets  Cr=[c.-c^x*^l.  effectively  cutting  off  the  objective 

value  of  (LP3)  at  0.  Ye  has  proved  that  this  version  of  the  algorithm  will 
converge,  but  is  polynomial  only  if  the  starting  point  w  is  "reasonably  close"  to 
the  optimal  solution.  Ye's  algorithm,  in  fact,  may  always  be  polynomial,  but 
this  seems  to  be  extremely  difficult  to  prove.  One  can  show  that  c^x*^*'  < 
c^x*^  for  the  cutting  objective  method,  which  implies  that  Cp  is  a  descent 

direction.  In  practice,  his  method  works  well  and  was  the  method  used  in  my 
implementation.  The  method  also  admits  an  elegant  convergence  criterion: 
(l!cpl|/|c^w|)<c.  For  my  implementation,  a  value  of  €=10"^  was  used.  Later  it 

will  be  shown  why  this  convergence  criterion  has  special  significance  in  some 
other  way. 

Given  a  linear  program  with  upper  bounds,  the  NPA  can  be  modified  to 
handle  them.  Assume  the  linear  program  is  in  the  form: 

min  c^x 

(LP5)  s.t.  Ax=b 

0  <  X  <  u 

if  lower  bounds  are  present  as  well,  the  above  form  can  be  obtained  by  a 
simple  translation  of  the  variables. 

Karmarkar  has  suggested  using  a  flipping  technique  to  handle  the  upper 
bounds.  As  the  algorithm  progresses,  if  x*|  >  (Uj/2).  then  Xj  is  flipped  by 

setting  x’|=Uj-x’|,  negating  the  j^  columns  of  c^  and  A,  and  adjusting  the  value 
of  b,  so  that  the  new  x  remains  feasible.  The  upper  bound  for  Xj  is  tested 
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using  the  ratio  test: 


«  =  max  \ 
j 


U  “X^ 
Uj  Xj 


(UjCp4,  -  xlISjXn+Dr 


provided  (UjCp,+  |  -  x*|cj)  ^  ®  r  • 


As  the  algorithm  converges,  the  "flipping"  of  the  variables  should  become  less 
frequent.  Note  that  this  procedure  preserves  the  upper  bounds  and  is  easy  to 
implement.  It  should  be  noted  that  a  proper  proof  of  convergence  has  not  yet 
been  exhibited.  The  procedure  has  worked  well  on  a  small  number  of  test 
problems. 
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Section  The  calculation  of  the  projected  gradient 


The  largest  computational  burden  In  any  algorithm  based  on  Karmarkar’s 
original  method  is  in  the  calculation  of  the  projected  gradient 
Cp=(l-B^(BB^r'B)DCp  during  every  iteration.  The  computation  Involved  can  be 

viewed  in  many  ways,  and.  in  particular,  as  the  following  least-squares 
problem: 

(LSI)  min  |DCr  -  B‘r(y)|2  . 

whose  solution  satisfies  the  normal  equations  (BB^)(J[)=BDCr.  It  follows  that 

Cp=DCp-B^(J[)  is  the  residual  vector  for  the  least-squares  problem  (LSI).  If 

(LP1)  has  a  non-degenerate  optimal  solution,  then  it  converges  to  the  optimal 
dual  solution.  In  any  case,  p  always  converges  to  zero.  In  the  case  of  the 
cutting  objective  method  of  Ye.  one  can  show  that  p=0  on  every  iteration. 
With  that  assumption,  one  can  then  write  Cp=D(lc,-c^x*^r-[A.-brir*^)  on  the 

iteration.  One  can  then  see  that  the  first  n  components  of  Cp  are  related  to 

the  complementary  slackness  conditions  that  hold  at  optimality  of  (LPI), 
because  of  the  construction  of  0.  Furthermore,  the  last  component  of  Cp 

corresponds  to  the  duality  gap  between  (LPI)  and  its  dual.  So  If  D  is 
converging  to  a  non-degenerate  solution,  the  size  of  the  norm  of  Cp  compared 

to  the  initial  estimate  (with  7t=0)  Is  a  reasonable  convergence  criterion.  In 
the  case  of  a  degenerate  optimal  solution,  the  convergence  criterion  is  still 
valid,  but  there  is  no  reason  to  believe  that  n  will  converge  to  a  specific  dual 
solution.  Since  a  degenerate  optimal  solution  corresponds  to  multiple  dual 
optima,  the  method  used  to  compute  Cp  will  determine  the  nature  of  the 

convergence  of  n.  However  the  norm  of  Cp  will  still  converge  to  zero.  ' 

M.  A.  Saunders  has  suggested  an  incremental  scheme  to  find  Cp  by  solving 

for  a  correction  term  to  the  vector  tt.  This  method  will  be  presented  in  the 
context  of  the  NPA  algorithm  presented  in  Section  2. 


Before  the  first  iteration,  we  set 
d  «-  Cp  and 

7t  ♦-  0  . 

On  each  iteration,  we  solve  the  least-squares  problem 
(LS2)  min  |Dd  -  bT(S)|j 

for  the  vector  (*). 

We  then  update  our  past  solutions: 

7t «-  n+S  , 
d  ^  d-[A,-bf  8  , 

and  calculate  the  final  projected  gradient: 

Cp^  Dd-pe^^+l  . 

Note  that  Dd=D(Cf.-lA,-brTr)  should  be  converging  to  zero,  by  the 

complementary  slackness  condition.  Hence,  we  are  solving  for  smaller 
correction  terms  to  tt  as  the  algorithm  proceeds. 

When  Ye's  cutting  objective  method  is  used,  it  is  also  necessary  to  update 
the  last  component  of  d  on  each  Iteration  before  solving  the  least-squares 
problem  (LS2). 


The  new  projective  algorithm  described  above  was  implemented  with  the 
cutting  objective  method  of  Ye  using  FORTRAN  on  an  IBM  3081  Model  KX2. 
Since  the  cutting  objective  method  was  used,  no  knowledge  of  the  final 
objective  value  was  needed  by  the  program  to  obtain  a  solution.  The  majority 
of  the  test  problems  were  the  same  as  those  run  by  Tomlin  [81  in  his  tests.  An 
additional  test  problem  was  created  from  an  exponential  counterexample  to  the 
simplex  method  and  came  from  Blair  Ill. 

In  order  to  calculate  Cp.  we  used  the  Iterative  algorithm  LSQR  of  Paige 

and  Saunders  [61.  The  incremental  method  described  in  the  previous  section 
was  used,  since  the  diminishing  vectors  Dd  in  the  least-squares  problem  (LS2) 
allow  LSQR  to  converge  more  rapidly  than  it  does  when  applied  to  (LSI).  A 
minimal  preconditioner  which  scaled  the  columns  of  each  least-squares 
problem  to  have  a  unit  Euclidean  norm  was  also  implemented.  Without  this, 
the  LSQR  algorithm  failed  to  converge  In  reasonable  time  on  some  problems. 

The  computational  results  are  presented  in  Table  1.  All  CPU  times 
represent  the  time  in  seconds  to  read  in  the  data  from  MPS  format,  solve  the 
linear  program,  and  write  out  a  solution.  Columns  labeled  MINOS  refer  to  the 
Iteration  counts  and  CPU  times  to  solve  the  same  problems  using  software 
based  on  the  Simplex  Method  developed  by  M.A.  Saunders  of  the  Systems 
Optimization  Laboratory  at  Stanford.  The  times  to  input  the  data  and  write 
out  the  solution  for  the  NPA  algorithm  and  for  MINOS  are  equal  because  MINOS 
subroutines  were  used  for  these  purposes. 

While  the  CPU  times  reflect  the  difficulty  of  using  LSQR  directly  in  the 
NPA  algorithm,  the  Iteration  counts  indicate  much  promise  for  NPA.  In  order 
to  reduce  the  work  done  on  each  iteration,  a  very  good  preconditioner  would 
have  to  be  Implemented.  It  should  be  noted  that  the  problems  BRANDY,  E226, 
and  BANOM  ail  have  null  variables  and  result  in  higher  iteration  counts  in  Phase 
I.  The  problem  ISRAEL  requires  a  large  amount  of  time  because  of  severe 
ill-conditioning.  This  problem  was  a  notable  omission  from  the  results 
presented  by  Karmarkar  in  Boston  [51. 

It  is  interesting  to  note  the  amount  of  work  done  by  this  rather  naive 
application  of  LSQR  to  find  Cp.  The  number  of  iterations  performed  by  LSQR 

can  be  linked  to  how  well-conditioned  the  least-squares  problems  are.  This 
conditioning  can  be  poor  on  either  the  first  few  iterations  of  Phase  I  or  the 


final  iterations  of  Phase  2.  depending  on  the  nature  of  the  problem.  Table  2 
summarizes  the  results  on  the  amount  of  work  performed  by  LSQR.  On  most 
problems,  the  number  of  iterations  performed  by  LSQR  on  (LS2)  seemed  to 
stabilize  as  the  projective  algorithm  converged.  This  usually  occurred  when 
the  algorithm  had  few  degeneracies  at  optimality. 

From  these  results,  it  is  apparent  that  success  or  failure  of  Karmarkar's 
algorithm  and  its  variants  will  depend  solely  on  how  fast  Cp  can  be  calculated. 

The  ill-conditioning  of  the  problems  may  cause  difficulties  at  the  beginning  of 
Phase  I  or  the  end  of  Phase  2.  Most  likely,  any  method  that  finds  Cp  will  find 

difficulties  because  of  this  ill-conditioning.  It  is  clear  that  further  research 
is  necessary  to  understand  how  Karmarkar’s  algorithm  creates  these 
conditioning  problems. 

It  is  also  important  to  note  the  low  iteration  counts  of  the  NPA 
algorithm,  when  the  cutting  objective  method  of  Ye  is  used.  The  promise 
provided  by  these  low  counts  is  encouraging  enough  to  warrant  further  research 
into  Karmarkar's  algorithm  and  its  variants. 
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Karmarkar's  algorithm  has  variants  which  are  easier  to  understand  in  the 
context  of  classical  linear  programming  theory.  These  variants  allow  one  to 
implement  computer  programs  that  run  with  low  Iteration  counts  and  solve 
linear  programs  with  upper  bounds. 

Karmarkar  claims  there  are  some  problems  which  will  be  solved  much 
faster  by  his  algorithm  than  by  the  standard  simplex  method,  but  it  is  probable 
that  an  adapatation  of  the  simplex  method  which  does  some  interior  probing 
may  also  do  well  on  these  same  problems.  For  most  linear  programs,  however, 
the  cost  of  doing  each  iteration  might  well  outweigh  any  low  iteration  counts 
obtainable.  It  is  still  questionable  whether  Karmarkar’s  algorithm  or  some 
variant  will  become  "the  method  of  choice"  for  linear  programming  in  the  near 
future. 

The  author  would  like  to  thank  George  Dantzig,  Narendra  Karmarkar, 
Edward  Klotz,  Nimrod  Megiddo,  and  Michael  Saunders  for  suggestions,  insights 
and  encouragement  as  this  work  progressed. 
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Foundation  Graduate  Fellowship.  Any  opinions,  findings,  conclusions  or 
recommendations  expressed  in  this  publication  are  those  of  the  author  and  do 
not  necessarily  reflect  the  views  of  the  National  Science  Foundation. 


Table  1=  Computational  Results  of  New  Projective  Algorithm. 
All  times  in  CPU  seconds  on  an  IBM  3081  Model  KX2. 


Problem 

Name 

Rows 

Columns 

MINOS 

(Simplex) 

New  Projective  Algorithm 

Total 

Itns. 

Phase  1 
Itns. 

Phase  II 
Itns. 

CPU 

Time 

EXPl 

10 

17 

52 

0.6 

6 

1 1 

0.4 

AFIRO 

28 

32 

6 

0.5 

3 

1 1 

0.8 

ADLITTLE 

57 

97 

126 

1.0 

9 

20 

12.3 

SHARE2B 

99 

79 

91 

I.O 

14 

67.4 

ISRAEL 

175 

142 

338 

4.2 

24 

636 

BRANDY 

221 

249 

292 

4.1 

19 

215 

E226 

226 

282 

572 

7.9 

■n 

42 

644 

BANDM 

306 

472 

362 

6.4 

u 

37 

Table  2-  LSQR  Iteration  Counts  for  NPA  algorithm. 

NS  =  Never  Stabilized,  continuing  to  increase. 


Problem 

Name 

LSQR  Iteration  Counts 

Max 

Attained 
Phase/ itn 

Minimum 

Maximum 

"Stable" 

EXPl 

12 

22 

12 

1/1 

AFIRO 

26 

55 

30 

1/1 

ADLITTLE 

83 

329 

~160 

1/1 

SHARE2B 

607 

786 

~775 

2/10 

ISRAEL 

270 

8853 

~650 

1/1 

BRANDY 

136 

1008 

NS 

2/18 

E226 

130 

1332 

-1100 

1/1 

BANDM 

175 

1679 

NS 
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